Calculate climate trends

Intro

This scripts calculates trends in biomass and climate variables from climate-agnostic sdm-predictions (“04-fit-sdms-random.qmd”)

Load packages & source functions

# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "RCurl", "devtools", "patchwork",
          "viridis", "RColorBrewer", "here", "raster", "concaveman") 

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library, character.only = T))

# Set path
home <- here::here()

# Source code for map plots
# For some reason I can't source it directly from github...
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
# devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
source(paste0(home, "/R/functions/map-plot.R"))
options(ggplot2.continuous.colour = "viridis")

Read data

# Read prediction grids
d <- read_csv(paste0(home, "/output/data_pred_grids_04_random_sdms.csv")) |>
  separate(group, c("species", "life_stage"), sep = "_", remove = FALSE)

Rows: 5660568 Columns: 42
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (6): substrate, month_year, ices_rect, group, species, life_stage
dbl (36): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Trim the prediction grid by species cumulative biomass

# Cumulative filter
# d_test <- data.frame(y = round(rnorm(10, 10, 2))) |> 
#   arrange(desc(y)) |>
#   mutate(cumsum = cumsum(y),
#          sum = sum(y),
#          perc = cumsum / sum)
# 
# d_test
# 
# d_test |> filter(perc < 0.95)

d_grid_keep <- d |>
  mutate(group = paste(group, quarter, sep = "_")) |> 
  group_by(X, Y, group) |> 
  # summarise across all years, but within group and quarter
  summarise(grid_mean = mean(est)) |>
  ungroup() |> 
  group_by(group) |>
  arrange(desc(grid_mean)) |> 
  mutate(cumsum = cumsum(grid_mean),
         sum = sum(grid_mean),
         perc = cumsum / sum) |> 
  ungroup() |> 
  filter(perc < 0.95) |> 
  mutate(id = paste(X, Y, group, sep = "_"))

mutate: changed 5,660,568 values (100%) of 'group' (0 new NA)

group_by: 3 grouping variables (X, Y, group)

summarise: now 195,192 rows and 4 columns, 2 group variables remaining (X, Y)

ungroup: no grouping variables

group_by: one grouping variable (group)

mutate (grouped): new variable 'cumsum' (double) with 195,192 unique values and 0% NA

                  new variable 'sum' (double) with 12 unique values and 0% NA

                  new variable 'perc' (double) with 195,181 unique values and 0% NA

ungroup: no grouping variables

filter: removed 21,973 rows (11%), 173,219 rows remaining

mutate: new variable 'id' (character) with 173,219 unique values and 0% NA

# Filter the main dataset to include these grid cells only
pred_grid_trim <- d |> 
  mutate(group = paste(group, quarter, sep = "_")) |> 
  mutate(id = paste(X, Y, group, sep = "_")) |> 
  filter(id %in% unique(d_grid_keep$id))

mutate: changed 5,660,568 values (100%) of 'group' (0 new NA)

mutate: new variable 'id' (character) with 195,192 unique values and 0% NA

filter: removed 637,217 rows (11%), 5,023,351 rows remaining

# Filter id's that are in the trends data frame before summarizing... because that is also trimmed by trends-quantiles. 
# So first trim by quantiles

Calculate trends and trim quantiles.

# This is if I use a spatial trend model
biom_yr_slopes <- pred_grid_trim |>
  mutate(id = paste(X, Y, group, sep = "_")) %>%
  split(.$id) |>
  purrr::map(~lm(exp(est) ~ year, data = .x)) |>
  purrr::map_df(broom::tidy, .id = 'id') |>
  filter(term == 'year') |>
  dplyr::select(id, estimate) |>
  rename(biom_slope = estimate) |> 
  separate(id, sep = "_", into = c("X", "Y", "species", "life_stage", "quarter"), remove = FALSE) |>
  mutate(X = as.numeric(X),
         Y = as.numeric(Y))

mutate: no changes

filter: removed 173,219 rows (50%), 173,219 rows remaining

rename: renamed one variable (biom_slope)

mutate: converted 'X' from character to double (0 new NA)

        converted 'Y' from character to double (0 new NA)

temp_yr_slopes <- pred_grid_trim |>
  mutate(id = paste(X, Y, group, sep = "_")) %>%
  split(.$id) |>
  purrr::map(~lm(temp ~ year, data = .x)) |>
  purrr::map_df(broom::tidy, .id = 'id') |>
  filter(term == 'year') |> 
  dplyr::select(id, estimate) |> 
  rename(temp_slope = estimate) |> 
  separate(id, sep = "_", into = c("X", "Y", "species", "life_stage", "quarter"), remove = FALSE) |> 
  mutate(X = as.numeric(X),
         Y = as.numeric(Y))

mutate: no changes

filter: removed 173,219 rows (50%), 173,219 rows remaining

rename: renamed one variable (temp_slope)

mutate: converted 'X' from character to double (0 new NA)

        converted 'Y' from character to double (0 new NA)

oxy_yr_slopes <- pred_grid_trim |>
  mutate(id = paste(X, Y, group, sep = "_")) %>%
  split(.$id) |>
  purrr::map(~lm(oxy ~ year, data = .x)) |>
  purrr::map_df(broom::tidy, .id = 'id') |>
  filter(term == 'year') |> 
  dplyr::select(id, estimate) |> 
  rename(oxy_slope = estimate) |> 
  separate(id, sep = "_", into = c("X", "Y", "species", "life_stage", "quarter"), remove = FALSE) |> 
  mutate(X = as.numeric(X),
         Y = as.numeric(Y))

mutate: no changes

filter: removed 173,219 rows (50%), 173,219 rows remaining

rename: renamed one variable (oxy_slope)

mutate: converted 'X' from character to double (0 new NA)

        converted 'Y' from character to double (0 new NA)

phi_yr_slopes <- pred_grid_trim |>
  mutate(id = paste(X, Y, group, sep = "_")) %>%
  split(.$id) |>
  purrr::map(~lm(phi ~ year, data = .x)) |>
  purrr::map_df(broom::tidy, .id = 'id') |>
  filter(term == 'year') |> 
  dplyr::select(id, estimate) |> 
  rename(phi_slope = estimate) |> 
  separate(id, sep = "_", into = c("X", "Y", "species", "life_stage", "quarter"), remove = FALSE) |> 
  mutate(X = as.numeric(X),
         Y = as.numeric(Y))

mutate: no changes

filter: removed 173,219 rows (50%), 173,219 rows remaining

rename: renamed one variable (phi_slope)

mutate: converted 'X' from character to double (0 new NA)

        converted 'Y' from character to double (0 new NA)

# Left join these data!
trends <- left_join(temp_yr_slopes, biom_yr_slopes, by = c("id", "X", "Y", "species", "life_stage", "quarter")) |> 
  left_join(oxy_yr_slopes, by = c("id", "X", "Y", "species", "life_stage", "quarter")) |> 
  left_join(phi_yr_slopes, by = c("id", "X", "Y", "species", "life_stage", "quarter"))

left_join: added one column (biom_slope)

           > rows only in x         0

           > rows only in y  (      0)

           > matched rows     173,219

           >                 =========

           > rows total       173,219

left_join: added one column (oxy_slope)

           > rows only in x         0

           > rows only in y  (      0)

           > matched rows     173,219

           >                 =========

           > rows total       173,219

left_join: added one column (phi_slope)

           > rows only in x         0

           > rows only in y  (      0)

           > matched rows     173,219

           >                 =========

           > rows total       173,219

# Trim slopes by quantiles
trends <- trends |> 
  group_by(species, life_stage, quarter) |> 
  mutate(upr_temp = quantile(temp_slope, probs = 0.995),
         lwr_temp = quantile(temp_slope, probs = 0.005),
         
         upr_oxy = quantile(oxy_slope, probs = 0.995),
         lwr_oxy = quantile(oxy_slope, probs = 0.005),
         
         upr_biom = quantile(biom_slope, probs = 0.995),
         lwr_biom = quantile(biom_slope, probs = 0.005),
         
         upr_phi = quantile(phi_slope, probs = 0.995),
         lwr_phi = quantile(phi_slope, probs = 0.005)) |>
  ungroup() |> 
  filter(temp_slope > lwr_temp & temp_slope < upr_temp) |>
  filter(oxy_slope > lwr_oxy & oxy_slope < upr_oxy) |> 
  filter(biom_slope > lwr_biom & biom_slope < upr_biom) |> 
  filter(phi_slope > lwr_phi & phi_slope < upr_phi)

group_by: 3 grouping variables (species, life_stage, quarter)

mutate (grouped): new variable 'upr_temp' (double) with 12 unique values and 0% NA

                  new variable 'lwr_temp' (double) with 12 unique values and 0% NA

                  new variable 'upr_oxy' (double) with 12 unique values and 0% NA

                  new variable 'lwr_oxy' (double) with 12 unique values and 0% NA

                  new variable 'upr_biom' (double) with 12 unique values and 0% NA

                  new variable 'lwr_biom' (double) with 12 unique values and 0% NA

                  new variable 'upr_phi' (double) with 12 unique values and 0% NA

                  new variable 'lwr_phi' (double) with 12 unique values and 0% NA

ungroup: no grouping variables

filter: removed 1,746 rows (1%), 171,473 rows remaining

filter: removed 1,746 rows (1%), 169,727 rows remaining

filter: removed 1,660 rows (1%), 168,067 rows remaining

filter: removed 245 rows (<1%), 167,822 rows remaining

Continue calculating average cell values

# filter the area-occupied-trimmed grid by the quantile-trimmed trends data frame
pred_grid_trim <- pred_grid_trim |> 
  filter(id %in% unique(trends$id))

filter: removed 156,513 rows (3%), 4,866,838 rows remaining

# ... Now that we have the final set of grid cells by species (after trimming by core area occupied and slope-quantiles), summarize grid cell average oxygen, temperature, and biomass
pred_grid_sum <- pred_grid_trim |> 
  group_by(id) |> 
  summarise(mean_log_biomass = mean(est),
            mean_temp = mean(temp),
            mean_oxy = mean(oxy),
            mean_phi = mean(phi)) |> 
  # Center the average across year variables
  mutate(mean_log_biomass_ct = mean_log_biomass - mean(mean_log_biomass),
         mean_temp_ct = mean_temp - mean(mean_temp),
         mean_oxy_ct = mean_oxy - mean(mean_oxy),
         mean_phi_ct = mean_phi - mean(mean_phi)) |>
  ungroup() |> 
  dplyr::select(id, mean_log_biomass, mean_temp, mean_oxy, mean_phi, 
                mean_log_biomass_ct, mean_temp_ct, mean_oxy_ct, mean_phi_ct)

group_by: one grouping variable (id)

summarise: now 167,822 rows and 5 columns, ungrouped

mutate: new variable 'mean_log_biomass_ct' (double) with 167,822 unique values and 0% NA

        new variable 'mean_temp_ct' (double) with 31,997 unique values and 0% NA

        new variable 'mean_oxy_ct' (double) with 31,996 unique values and 0% NA

        new variable 'mean_phi_ct' (double) with 167,822 unique values and 0% NA

ungroup: no grouping variables

# Join in the average variables to the trends data
trends <- trends |> 
  mutate(id = paste(X, Y, species, life_stage, quarter, sep = "_")) |> 
  left_join(pred_grid_sum, by = "id")

mutate: no changes

left_join: added 8 columns (mean_log_biomass, mean_temp, mean_oxy, mean_phi, mean_log_biomass_ct, …)

           > rows only in x         0

           > rows only in y  (      0)

           > matched rows     167,822

           >                 =========

           > rows total       167,822

# Scale variables: which group?
trends <- trends |>
  mutate(group = paste(species, life_stage, sep = "_")) |>
  group_by(group) |>
  mutate(temp_slope_sc = (temp_slope - mean(temp_slope)) / sd(temp_slope),
         oxy_slope_sc = (oxy_slope - mean(oxy_slope)) / sd(oxy_slope),
         phi_slope_sc = (phi_slope - mean(phi_slope)) / sd(phi_slope)) |>
  ungroup() |>
  mutate(quarter_f = as.factor(quarter))

mutate: new variable 'group' (character) with 6 unique values and 0% NA

group_by: one grouping variable (group)

mutate (grouped): new variable 'temp_slope_sc' (double) with 167,822 unique values and 0% NA

                  new variable 'oxy_slope_sc' (double) with 167,822 unique values and 0% NA

                  new variable 'phi_slope_sc' (double) with 167,822 unique values and 0% NA

ungroup: no grouping variables

mutate: new variable 'quarter_f' (factor) with 2 unique values and 0% NA

# Fix names
trends <- trends |>
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage)) |>
  mutate(group = paste(species, life_stage, sep = "_"))

mutate: changed 167,822 values (100%) of 'species' (0 new NA)

        changed 167,822 values (100%) of 'life_stage' (0 new NA)

mutate: changed 167,822 values (100%) of 'group' (0 new NA)

Plot trends!

# Because we trim the grids, let's add a polygon that is the overall outside border.
border_df <- trends
z <- concaveman(cbind(border_df$X, border_df$Y), concavity = 2)
plot(z)

z_poly <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(z)), ID = 1)))
domain <- as.data.frame(z_poly@polygons[[1]]@Polygons[[1]]@coords)
colnames(domain) <- c("X", "Y")

ggplot(domain, aes(x = X, y = Y)) +
  geom_polygon(fill = NA, color = "tomato")

# Plot velocities!
long_trends <- trends |> 
  # Not using phi here!
  pivot_longer(c("biom_slope", "temp_slope", "oxy_slope")) |> 
  #pivot_longer(c("biom_slope", "temp_slope", "oxy_slope", "phi_slope")) |> 
  mutate(name2 = ifelse(name == "biom_slope", "Biotic", NA),
         name2 = ifelse(name == "temp_slope", "Temperature", name2),
         name2 = ifelse(name == "oxy_slope", "Oxygen", name2)#,
         #name2 = ifelse(name == "phi_slope", "Metabolic index", name2)
         ) |> 
  mutate(loop_id = paste(paste0("Q", quarter), name2, sep = " ")) |> 
  mutate(unit = "(kg/decade)",
         unit = ifelse(name2 == "Temperature", "(°C/decade)", unit),
         unit = ifelse(name2 == "Oxygen", "(ml/L/decade)", unit)#,
         #unit = ifelse(name2 == "Metabolic index", "(M.I./decade)", unit)
         )

ggplot(long_trends, aes(value)) + 
  facet_wrap(~name, scales = "free") +
  geom_histogram()

# Loop through the variables... Note if oxygen, the tail is negative so ... trim by a small quantile

for(i in unique(long_trends$loop_id)) {
  
  dd <- filter(long_trends, loop_id == i)
  
  print(
    plot_map_fc + 
      geom_raster(data = dd, aes(X*1000, Y*1000, fill = value)) +
      facet_grid(life_stage ~ species) +
      scale_fill_viridis(option = "rocket", name = paste(i, "trend", dd$unit[1])) +
      geom_polygon(data = domain, aes(X*1000, Y*1000),
                   fill = NA, color = "red", size = 0.2)
    )
    
  # ggsave(paste0(home, "/figures/supp/trend_", str_replace(i, " ", "_"), ".pdf"), width = 17, height = 13, units = "cm")

}

Save data

# Save data
write_csv(trends, paste0(home, "/data/clean/trends.csv"))